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ABSTRACT 



The problem of narrow band interference while transmitting broad-band signals 
like Direct Sequence Spread Spectrum is a common source of problems in Electronic 
Warfare. This can occur either due to intentional jamming or due to unavoidable 
signal sources present in the vicinity of the receiver. Lack of improper information 
on these narrow band interferers makes it difficult to cancel them. 

In this thesis the above problem is addressed by using an adaptive notch filtering 
technique. Before adopting such a technique other methods like the Least Square 
Estimator and the Maximum Likelihood Estimator were explored. However the Kwan 
and Martin adaptive notch filter structure was found both relevant and suitable for 
the problem of interest. The Kwan and Martin method has the difficulty of increasing 
hardware complexity with number of notches. This makes it difficult to implement 
in real time. A new algorithm was developed for the purpose of implementing the 
structure in real-time. This new algorithm offers the same performance at reduced 
hardware complexity. This algorithm was simulated and the results were presented. A 
hardware feasibility is discussed by proposing a simple structure based upon existing 
commercial signal processing chips. 
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I. INTRODUCTION 



A. MOTIVATION 

The problem of suppressing narrow-band noise from a broad-band spectrum is 
of considerable importance in the areas of Electronic Warfare and Anti Submarine 
Warfare scenarios. In this thesis a workable solution is proposed for the problem cre- 
ated by narrow-band interference. This solution will come in the form of an adaptive 
digital notch filter. However before the details of the solution to the problem are 
discussed, a brief description of the problem is given 

1. Electronic Warfare(EW) 

A typical EW scenario is given in Fig 1.1. In this scenario a transmitter 
transmits information over a long distance. As an Electronic Counter Measure(ECM) 
this information is coded and transmitted as a Direct Sequence Spread Spectrum 
(DSSS) signal. Inherently this signal is a broad-band signal. However at the EW 
receiver there is often narrow-band interference from such things as push to talk 
systems(PTS) that swamp out the received DSSS signal. Sometimes for reasons of 
signal security PTS operating frequencies are varying with time. Under conditions 
such as these the EW receiver cannot function effectively. For proper functioning of 
the EW receiver we must enhance the received signal by selectively and adaptively 
suppressing these narrow-band signals. 

2. Assumptions 

This thesis addresses the problem arising at the tactical data communica- 
tion link. In this EW scenario as dipicted in Fig 1.1 we are attempting to perform 
signal analysis on a DSSS signal using a EW reciever. This EW receiver is not the des- 
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ignated receiver and hence does not have the spreading code available. For this reason 
the narrow-band interference severely limits the ability to determine the charaterstics 
of the recieved signal such as carrier frequency, chip-frequency etc. 

A typical DSSS system is the Tactical Information Distribution Systems 
(JTIDS)[Ref. 7]. Data on such a link is generally in the form of digital messages. 
A typical L-band (960-1215 Mhz) JTIDS allows transfer of digital data between any 
properly equipped users within line-of-sight. The typical RF band-width required is 
around lOMhz. 

Assuming a Superheterodyne EW-receiver, the band-width of the Interme- 
diate Frequency (IF) need be only lOMhz. This enables us to digitize the signal at 
the base-band with a (20-50)Mhz Analog to Digital(A/D) converter. It is assumed 
that this signal is digitized and converted into a floating point number. In this thesis 
it is assumed that the floating point arithmetic is of sufficient precision that rounding 
errors are not a significant problem. 

B. PROBLEM DEFINITION 

The signal described in section (I. A. 3) can be thought of as a broad-band signal 
with additive narrow’-band signals occurring at arbitrary frequencies and times. These 
narrow-band signals can be modeled as either pure sinusoidal signals, as outputs of a 
narrow-band filter excited by white noise, or as any band-limited narrow’-band signal. 

We shall assume that the signal yk is received by the receiver and is defined as 

3 

Vk = m + x\ + 7* (1.1) 

<=] 

where the interferer, is either a pure sinusoid or a narrow-band signal whose 
frequency (i.e center frequency) can vary slowly with time. The term Wk is the 
broad-band DSSS signal of interest and 7 * is White Gaussian Noise(WGN) with zero 
mean and variance a 2 ,/V(0, a 2 ) The goal is to remove only the narrow-band interferes 
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and obtain the output (w^ + 7*) for further processing. A typical spectrum for y * is 
given in Fig 1.2. 

In this thesis, the design for an adaptive digital filter that will achieve the goal 
of removing narrow-band interferes from a broad-band signal of interest is developed. 
This thesis is organized into five chapters, chapter I is the introduction. The chapter 
II concentrates on building necessary background on such things as the Least Square 
Estimator and Maximum- Likelihood Estimator. Chapter III briefly explains existing 
methods of adaptive filtering with emphasis on adaptive notch filtering techniques 
. The new algorithm [Ref. 26] developed for the purpose of addressing the current 
problem of interest is also explained. Chapter IV is totally devoted to modelling 
various signals required for testing. Simulation results are discussed and presented. 
Chapter V provides a brief survey of existing hardware components and a look at 
possible hardware structures and their limitations. Chapter VI is conclusions. 
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Fig 



II. ESTIMATION TECHNIQUES 



A. DIGITAL FILTERS 

Digital Filters can be characterized by their impulse response, their trans- 
fer function, or by a difference equation [Ref. 13] . A typical Infinite-Impulse- 
Response(IIR) filter represented in difference equation form is given by 

x k = aix k -i + a 2 x k . 2 + a 3 x k . 3 + + 6 2 «fc-i (2.1) 

while an Finite Impulse Response(FIR) filter is given by 

x k = biU k + b 2 u k - 1 + b 3 u k - 3 (2.2) 

where u k is the input to the filter and x k is the output of the filter. x k and u k are 
the sampled values of the continuous functions x(t) and u(t). By choosing T, as the 
sampling time we notionally write the discrete signal as 

x k = x(k) = x(kT a ) (2.3) 

Let A'(ff) be the spectrum of the signal x k , where is the angular frequency 
component present in the actual signal. For analysis purposes consider normalized 
frequency given a s« = f IT,. With this definition we do not have to use the actual fre- 
quencies but only the normalized frequencies. Assuming that there is no aliasing while 
sampling the maximum normalized frequency is 0.5 cycle/sample or n rad/sample. 

1. Some Hardware Schemes 

Fig 2.1 represents a conventional Direct Form II representation of a digital 
filter [Ref. 18]. However for high-speed throughput it is often necessary to have a 
pipelined architecture. FIR filters are highly amenable for pipelining thus giving a 
high throughput rate for computations. 
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Fig 2.1 



2. Modelling Hardware Multiplier and Adder 

Multiplication of two variables say bj and u k in hardware can be modelled 
as 

b\*Uk = bi * UkDj = b\UkD\ (2.4) 

where * is the hardware multiplier and * is the ideal multiplier and D\ is the propa- 
gation delay associated with the multiplication. Similarly a hardware adder is desig- 
nated as -f while the ideal adder is given as + and the corresponding delay associated 
with the adder is D 2 . Then an addition of two variables u k and u k ~j can be written 
cis 

Uk+Uk-i = (ujt + u k -\)D 2 (2.5) 

Each of the previous devices can be incorporated into a pipeline digital filter structure 
by following it with a synchronizing register. The clock period must then be less than 
D + t s + t r where t, is the register setup time and t r the register propagation delay. 

3. Pipelining an FIR filter 

In this section we wish to realize 

x k = bi * u k + b 2 * u k - 1 + 63 * u k - 2 (2.6) 

It is acceptable to have an overall delay D‘ % , However it is desirable to minimize N: 

x k = {&] *u k + b 2 * Uk-\ + b 3 * u k - 2 } (2.7) 

First associate one D with each multiplier and form real multipliers by the notation 

D* = *. Then we get 

ik = T» A_1 {b]D * Uk + b 2 D * uk — 1 -f b 3 D * u k - 2 } (2.8) 

Now form the real multiplier 

Xk = {bi*Uk + b 2 *uk — 1 + b 3 *Uk- 2 } (2.9) 
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Now we use one additional D to associate with the adding of b 2 *Uk~i and b 3 *Uk- 2 . 
(Note: It is important to choose the most delayed values first in order to minimize 

A r ) 

ik — D** 2 {Db\iuk + D(b 2 *u k -, + 63 * A: — 2 ) } (2.10) 

Now form the real adder 

ik = D A 2 ^Dbiiuk + ^iuk-i + b^iuk-?)} (2-11) 

Now we use an additional delay D to associate with the final adder to get a real 
adder: 

ik = D N ~ 3 {Dhiuk+ihiuk-r+bsiuk-t)} (2.12) 

Finally we choose the system delay such that D = z~ l and replace u ^_ 2 = 
z~ 2 Uk and Uk~ 1 = z~ 1 Uk, thus we get 

Xk = z _A+3 (2.13) 

Now factor z -1 out of the expression to obtain 

i k = z~ N+2 {b 1 *Uk+(b 2 *Uk+z~ 1 b 3 *Uk)} (2.14) 



Note that iV = 2 is the minimum possible value of N for a causal filter. Choose 
N = 2 then 



i k = {bi*Uk+(b 2 *Uk+z l b 3 *Uk)} 



(2.15) 



Equation (2.15) is very convenient for pipelining and is given one to one in Fig 2.2. 

The above procedure is general for any FIR filter. Pipelining an FIR filter 
results in maximizing the sampling or throughput frequency at the expense of a delay 
or latency in the availability of the output. 
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4. Pipelining HR filters 

Pipelining an HR filter is considerably more difficult than an FIR filter. 
This is mainly because the delay in the availability of the output makes it impossible 
to feed back output values with short delay as is required for the straightforward HR 
implementation. As an example consider a simple 1st order HR filter 



x k = a l * x k .i + u k (2.16) 

As in the FIR case we introduce a delay D A as follows. 

x k = J9 A (ai * x k -i + u k ) (217) 

Now we use one delay to form the real multiplier: 

x k = T> A_1 (aiJ9 * x k -i + Du k ) (2.18) 

= D N ~ 1 {a i *Xk-i + Du k ) (2.19) 

We use one more D to form the real adder 

x k = D N ~ 2 (a l D * x k ^i+u k ) (2.20) 

Finally assume D = z~ l and x k ^i = z~ l x k then 

Xk = z-^faiz-'xk+z-'uk) ( 2 . 21 ) 

= z- N+ \ ai ix k +Uk) ( 2 . 22 ) 



For x k = x k , N must be zero. However > 1 is required for a causal filter. 

Solutions to this problem exist. They use a higher order difference equation 
representation of the filter with equivalent characteristics, so that feedback loop delay 
greater than 1 can be tolerated. More details can be found in the literature [Ref. 14] 
regarding pipelining of HR digital filters. 
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B. SIGNAL CHARACTERIZATION 



The signal yk of equation (1.1) can be characterized either in the time domain 
or frequency domain and the signal can be easily transformed from one domain to the 
other via a linear transformation. 

1. Time Domain Representation 

In the time domain representation, the signal y * is modelled as the output 
of a linear system which may be either an all pole or a pole-zero system [Ref. 13]. 
The input is assumed to be white noise but this input to the system is not accessible 
and is only conceptual in nature. Let the system under consideration be 

V 9 

Vh = a iyk-i *F (2.23) 

t=l j~0 

where 7 * is a white noise and yk is the output of the system. Equation (2.23) also 
can be represented as a Transfer Function(TF) [Ref. 17] 

"M = Wi (2 ' 24) 



or 




where 

M z ) = 1 - J2 a ' z ~' 

;=i 

and 

aw = EW’ 

j=0 

we define the parameter vector as 



(2.25) 



(2.26) 



(2.27) 



p ( — [aij <*2, • • • t a p > &2> • • • , M 



(2.28) 



The parameter p completely characterizes the signal x* The above system defined 
by the equation (2.23) is also known as an Auto-Regressive-Moving-Average(ARMA) 
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model [Ref. 12], [Ref. 15]. It is conventional to denote a p poles and q zeros system 
as ARMA(p.q). 

2. Frequency Domain representation 

The signal j/t of equation (2.23) can also be represented as a vector 

Y N = [yi---yN] (2.29) 



which can be transformed into the frequency domain by the DFT relation 




(2.30) 



where \V = e and j = In general gk is a complex quantity. However we 

restrict our interest to the power spectrum a real quantity given as 



s k = gk9k 



(2.31) 



The series 5 a t 



Sn = [si • • • sa'] 



(2.32) 



is another form of representing the signal and is designated as the power spectrum of 
the signal Xk- Even though the signal in this domain is very convenient to handle, 
computational complexity and other considerations limit its use for online application. 
Also in this representation phase information of the signal is lost. 



C. SPECTRUM ESTIMATION 
1. Least Square Estimator 

Consider the signal y k the output of a linear system excited by a white 
noise as given by the equation (2.23). Now our problem is to estimate the parameter 
vector p by obtaining successive measurements of yk- We shall define the vector 

= [j/t-i) ' ' ' i Pk-p> 7 k- > 7*-?] (2.33) 
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and Zk = yk so that the difference equation (2.23) can be rewritten as 

Zk = P fcXfc 

an elegant solution [Ref. 2] for the above problem can be written as 

P k = P*-i + 



where 



and 



-l 



H * = [Y x « x «- 



Ck = Zk — Pk-l*k 



(2.34) 



(2.35) 



(2.36) 



(2.37) 



(2.38) 



The block solution for the equation (2.34) can be also written in the form 

P* = H/. (£ XkZk'j 

Equation (2.37) can also viewed as Axk — £ 7 * and this is represented in Fig 2.3. 

Computing the matrix Hjt via equation (2.36) is computationally inefficient 
and a recursive solution via the matrix inversion lema [Ref. 3] is preferred and is given 
by 



Hjt = Hjt_i — 



(2.39) 



1 +xj.H jt _iXj t 

Staticians refer to the H^- matrix as the covariance matrix while optimization people 
refer to as the Hessian of the objective function, where the objective function is 
given by 

Jk = Ye] (2.40) 

i=i 

Appendix A gives a computer program based on this approach that was used for this 
thesis work to identify system parameters via Equations (2.35), (2.36), (2.37), (2.39). 
The above method is known to work well as long as the measured value of the state 
of the filter is free from noise. In the event Zk is not noise free, estimates are known 
to be biased. 
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Fig. 2.3 



2. Maximum Likelihood(ML) Estimator 

The ML estimator for the above problem is conventionally obtained by 
considering a log-likelihood function [Ref. 16] and minimizing it. However for getting 
a better physical understanding the approach given by Voting Ref. 2] will be followed. 
Consider a system as 



Vk = 



' BjzY 

.M z ). 



ik 



and an auxiliary system referred to as the model is given as 



(2.41) 



&k 



B(z) 



MW 



7k 



(2.42) 



and we define the function to be minimized a s L = e£( p) where e* = yk — Xk- 
To minimize this function the derivatives and the Hessian of the function L where 
L = (yk — £k) 2 need to be obtained. First take partial derivative of L with respect to 
each parameter. For example 



dL -0( * \( dik \- 0 ( dik \ 
dbi --{yk **)( db )~ 2ek ^ dh ) 



(2.43) 



The partial ||^ can be obtained by differentiating equation(2.42) as 



- 1 



dbi (A(^) 



6j b 2 z 



-1 



1 



lA(z) 



1 — aiz~ l — a 2 z ~ 2 — a^z~ z 
7 k 



Ik 



Similarly the other partial can be obtained as 

JM|W 1 = 

da, \ 



fl(-) -■ 

,[>4Wf 



7k 



Ik 



(2.44) 

(2.45) 



(2.46) 

(2.47) 
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since r^-{i?(i)} = 0 as B(z) has no terms containing aj. Using equations (2.45) and 

(2.47) the partials are: 



dik 


1 




db x Uk ~ 


M*). 


ik 


dx k 


B(z) 1 


da ! Xk ~ 




dy k 


1 




d ai Vk 


iMz)\ 


Vk 



by defining 

p! = [ Q i» ° 2 , Qa, 6 2 ] 

x l = [ x fc-l> x It-2i x I-3> u I> u i-l] 

Z 1 = [j/jt-n yk-2i Vk-31 U k> u fc-l) 

the Hessian may be approximated by: 

h ‘ = (!H 

Using the above equations the parameter can be estimated via 



(2.48) 

(2.49) 

(2.50) 

(2.51) 



(2.52) 

(2.53) 

(2.54) 



(2.55) 



P* = P*-i + (2.56) 

where /r is the step size in incrementing the parameter vector. A block schematic for 
the algorithm was given in Fig 2.4. A serious drawback of this method is the stability 
while incrementing the parameters. 

3. Other Methods 

There are many other methods such as the Instrumental Variable (IV) 
approach, the Approximate-Maximum-Likelihood (AML) Method and a combination 
method IV-AML available in the literature [Ref. 2]. 
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Maximum Likelyhood Method 
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Fig. 2.4 



D. SELECTION OF ESTIMATION TECHNIQUE 



Before we select the appropriate method for the problem of interest, model-order 
problems must be considered. 

1. Model Order Problems 

In the Least-Square Estimator for the parameter vector p of the signal 
Xk we need to assume the dimension of p or the order of the system. Since the 
system order is often not known a higher-order model estimate is assumed. Now we 
investigate higher-order model fits for a lower-order data. 

Consider the system equation (2.23) with parameter given as 

p' = [1.7,-1.53,0.648,1,0.6] (2.57) 

For the above system a Random Binary Noise(RBN) was given as the input Uk and 
the corresponding output Xk is obtained. For the time series Uk and Xk an ARMA(6,1) 
model was fit using the program given in the appendix and the parameter was esti- 
mated as 

p‘ = [2.26533,-2.79577,2.17249,-1.09877,0.446598,-0.110151,0.998175] (2.58) 

Fig 2.5 gives the parametric spectrum for the above. 

For the same data an ARMA(4,1) model was used and generating a pa- 
rameter vector 



p‘ = [2.16988,-2.41286,1.47436,-0.365412,0.999317] (2.59) 

Fig 2.6 gives the parametric spectrum for the above values. 

Lastly for the same data an ARMA(8,1) model was used and generating a 
parameter vector 

p ( = [2.276,-2.835,2.267,-1.282,0.703,-0.353,0.143,-0.0035,0.9999] (2.60) 
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Fig 2.7 gives the corresponding spectrum. 



It is interesting to compare these with the actual parametric spectrum given 
in Fig 2.8 For the sake of completion and comparison, the spectrum given in Fig 2.9 
of the output time series Xk computed via a DFT program is also included. This 
simulation demonstrates that the existence of multiple solutions, hence an important 
step in the estimation procedure is to choose an appropriate model order. 

2. Choice of the Method 

Since the filter hardware must be implemented in real-time at the frequen- 
cies of 10 to 20 Mhz computational complexities must be kept to a minimum. These 
methods although providing good performance, are not well suited to hardware im- 
plementations. This motivates looking into new filtering methods. 
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III. METHODS for FILTERING 



A. FILTERING TECHNIQUES 

Filtering in a broad sense is selectively suppressing a portion of the spectrum 
of the given signal. In this chapter several filtering techniques are explored in the 
light of the narrow-band interference problem in order to identify applicable filtering 
techniques. 

1. Filtering via FFT 

Any given signal can be conveniently transformed into its power spectrum 
via equation 2.31. Assuming the spectrum of the desired filter to be defined by w^. 
The weighted output is given by: 

Vk = SkWk (3.1) 

and the corresponding filtered output t/jt is obtained by taking the inverse DFT of the 
signal Vk . A simple block schematic representing this idea is given in Fig 3.1. Details 
of this approach are available in various references [Ref. 22]. 

Fast algorithms such as the FFT for computing the DFT make it possible 
to do the above process in real time by using dedicated hardware. Honeywell makes 
the HDSP66110 and HDSP66210 Digital Signal Processing chip pair w’hich are ideally 
suited for these applications. This chip pair can perform a single complex multiply in 
40ns [Ref. 28]. However the problem of filtering adaptively [Ref. 10] (i.e varying Wk 
according to a criterion) demands more computing power which can be obtained by 
augmenting the DSP chip pair with a processor like the R3000 which has a Reduced 
Instruction Set Computer(RISC) architecture with a high instruction execution rate 
of approximately 25 Million Instructions Per Second(MIPS) [Ref. 27]. 
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Fig. 3.1 



2. Recursive DFT 



Implementation of a higher order FIH filter using a recursive DFT [Ref. 
23], [Ref. 24] is also very convenient for eliminating narrowband interference. The 
basic concept is that tin* FIR filter is expressed as a product of two filter sections. 
One section is a filter with its zeros being equally spaced on the unit circle. 'I his 
is achieved by a delay. The second section is a pole-producing section. Pole zero 
cancellation results in the desired FIR filter. More details can be seen in t lie references 
[Ref. 23], [Ref. 24]. 

3. Adaptive Filtering 

Adaptive filters can be placed into four classes based upon the choice of 
the training sequence and the reference model for adaptation. Simon [Ref. 21]has 
classified the Adaptive Filters into the following four classes 

• Identification - Class 1 

• Inverse Modelling - Class 11 

• Prediction - Class 111 

• Interference Canceling - Class IV 

Fig [3.2 - 3.5] give the block schematics of the various classes of adaptive 
filters. In our problem, the reference rnodcl is naturally a narrow-band bandpass filter 
since the interference signal is a narrow-band signal. The class of filtering that best 
suits our problem is a combination of Class III and Class IV type of filter [I ig 3.4 
and 3.5]. This logically points to an adaptive notch filter. Due to apriori knowledge 
of the interference signal, only the parametric approach was. However the 1)1 1 
based techniques may also ofTer a good solution and should be the subject of future 
investigations. 
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B. ADAPTIVE NOTCH FILTERS 



Notch filters for removing multiple narrow-band interference can be categorized 
into four broad categories illustrated in Fig 3.6 through 3.9. The first two categories, 
Figures 3.6 and 3.7, are cascaded second order notches with each second-order section 
removing one frequency. The next two categories, 

Figures 3.8 and 3.9, are higher-order notches that eliminate multiple frequen- 
cies. In all of the categories, it is possible to use FIR filters (i.e all zero filters) which 
are easily pipelined and can be made truly linear phase. However, HR notch filters 
require substantially fewer multipliers and adders than FIR notch filters. Thus IIR 
pipelining may become an important issue. In this thesis we limit our discussion only 
to the first two categories illustrated in Figures 3.6 and 3.7. 

1. Second-Order Cascaded Notch Filters 

The second-order notch filter is used in cascade and in-line with the signal 
as shown in Figure 3.6. The transfer function for such a notch filter is given by Kwan 
and Martin [Ref. 8] as: 



H x (z) = 1 -H bp (z) 

[h (1 + *-*)(!-*-») 

\ 2 1 - (2 - k 2 - k\)z - 1 + (1 - k 2 )z~* 



2-fe 2 1 ~ 2(2 



(3.2) 

(3.3) 

(3.4) 



2 1 - (2 - k 2 - k\)z - 1 + (1 - k 2 )z ~ 2 

For arbitrary values of k\ and k 2 , this is a symmetric notch filter with unity 
gain at DC and the Nyquist frequency. If k 2 is kept constant, then the 3db notch 
width is also kept constant. Thus k\ may be adapted to remove one narrow-band 
signal. A cascade of such filters can be used to remove multiple narrow- band signals 
[Ref. 8]. Constants k\ and k 2 are related to the pole radius r and the normalized pole 
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Fig 3.9 



frequency O p as 



h = (l-r 2 ) (3.5) 

ki = y^l -f r 2 — 2rcos6~ p (3.6) 

It is important to bear in mind that the above transfer function has unity gain at 

Opeak = 2 arcsin < — > (3.7) 

which is different from 0 P [Ref. 8] 

2. Second-Order Cascaded Signal Canceler 

The cascaded second-order signal canceler approach shown in Figure 3.7 
has the advantage that the desired signal does not pass through the adaptive filter. 
Instead, the band-pass filter is used to detect the narrow-band signal which is then 
subtracted from the desired signal. A constant 3db bandwidth notch can be achieved 
by selecting a band-pass filter with the transfer function: 

%(.-) = 



Mi -»~ 2 ) 

2 D(z) 

A’(z) 

D(z) 



(3.8) 

(3.9) 



where D(z) is the same denominator as H^{z) in equation (3.4) The signal-canceler 
structure is also nice for adaptation because it is relatively easy to generate sensitivity 
functions which are related to the gradient of fhp(z) with respect to the frequency 
parameter fcj. Figure 3.10 shows the block diagram of an adaptive version of this 
filter. The sensitivity function s(n ) is obtained by differentiating the error signal 
e(n) with respect to the parameter k\. This can be easily derived by noting that 
e(?j) = x(n) — y(n) and 



de(n) dy(n ) 

dk ] = dJ~ 



(3.10) 
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to get the above partial derivative we consider 

y(n) = H bp (z)x(n) 
N(z) 



D(z) 



x(n) 



(3-11) 

(3.12) 



the above equation (3.12) can be easily differentiated with respect to ki by using the 
equation (2.47) to obtain: 



dy{n) 



dk, 

DW) 

dyjn) 
dh 



N(z) 



2 k x z 



-l 



x(n) 



(3,13) 



.PWf 

by recognizing that Hb p (z) = from equation (3.12) the above equation can be 
written as 



«bU) 



2kiz 



-l 



and by defining 



H„U) = 



DU) J 

2k x z~ l 

DU) 



x(n) 



(3.14) 



(3.15) 



we get the sensitivity function as 



s(n) = luT = = H bU)u„U) 



(3-16) 



The equation (3.16) is given as a block schematic in Fig 3.10. The parameter ki may 
then be adapted by the formula: 



k } (n + 1) = ki(ri) — fie(n)- 



4”) 

|s(n)|| 2 



(3.17) 



Computing ||s(n)|| 2 

Ideally ||s(n)|| 2 can be calculated by the relation 



Nn )|| 2 = S ( ? T 

i=n-N 



(3.18) 
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The above scheme computes the average of the past N samples. However a weighted 
average of the past samples with highest weight for the recent sample can be done 
quite easily by a first-order low pass filtering [s(n)] 2 as follows: 

v n = u„_iA + (1 - A)s 2 (n) (3.19) 



However the above lowpass filtering can also be carried out by a second order filter 



such as: 






6j(l — 2z 1 + z 2 ) 
1 — 2rjkz~ l + z ~ 2 



[s(n)] 2 



Vn — 



&i(l — 2cos(56)z 1 + 2 2 ) 
1 — 2rjkz~ l + z~ 2 



[s(n)] 2 



Yet another way to estimate ||.s(n)|| 2 i simply: 



(3.20) 

(3.21) 



V„ 



£ [*■'(«)]’ 



(3.22) 



the above form is defined as zero order forgetting. It should be noted that equation 
(3.21) was used by Kwan and Martin [Ref. 8]. 



C. MULTIPLE NARROWBAND SIGNAL SUPPRESSION 

The second-order implementations of section (1.) and (2.) offer considerable 
advantage both in hardware complexity and in adaptive performance .for multiple 
narrow band signal suppression [Ref. 8], [Ref. 25], [Ref. 26]. 

1. The Kwan and Martin Filter 

In a recent paper by Kwan and Martin [Ref. 8], the problem of detecting 
and enhancing sinusoidal signals in the presence of noise is addressed with a cascade 
of HR adaptive notch filters which are used to eliminate the sinusoids. Each of the 
sinusoids is eliminated by a bandpass filter whose output is an enhanced version of 
one of the sinusoids. Hence this remarkable structure can perform both tasks with a 
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single adaptive filter configuration which is shown to be highly robust and performs 
extremely well. 

The major disadvantage of the Kwan and Martin structure is that the 
number of biquad sections needed in the adaptive filter configuration is given by 
N(N+3)/2, where N is the number of sinusoids to be detected and removed. This 
becomes impractical in real-time situations with more than 4 sinusoids due to the 
geometric increase in the required hardware. 

a. Kwan and Martin Structure 

The Kwan and Martin structure consists of a cascade of HR notch 
filters one stage of which is shown in Figure 3.11. Each stage consists of a bandpass 
filter with zeros at DC and the Nyquist frequency and unity gain at its peak frequency 
u Such a filter w'ould have the following z-domain transfer function 




1 — r? 1 — z~ 2 

2 1 — 2r i cos0iZ~ 1 + r 2 z~ 2 



(3.23) 



where 



r, = pole radius of the i-th section 

6, = 2iriVi/u! s 

u>, = peak frequency of the i-th section 
u> 3 = sampling frequency 

Kw T an and Martin identify two different methods for adapting the 
filter. Most of their derivation is based on what they call the constant bandwidth filter 
in which the pole radius r, is a constant and only the frequency u>, is adapted. An 
alternative approach which keeps a constant Q is also discussed in Kwan and Martin 
paper. In addition, Kwan and Martin select the adaptive quantity in such a w r ay that 
it is fairly easy to determine the notch frequency from the adaptive parameter. From 
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Figure 3.10, we see that the notch filter for each section is the difference between 1 
and the bandpass filter, hence: 



Hh(z) = 1 - H‘ h (z) (3.24) 

b. Calculation of the gradient 

The basic structure of the Kwan and Martin adaptive filter shown in 
Figure 3.12 is a cascade of N sections of the form of Figure 3.7. The overall transfer 
function is given by 



T(z) = rW( 2 ) (3-25) 

1=1 

= n (i - k ) < 3 - 26 ) 

t=l 

Kwan and Martin choose as their objective function J(z ) the square 
of the output of the final stage of the cascade: 



J(-) = WWf 

= [rwflA'M] 2 



(3.27) 

(3.28) 



Hence, the gradient of the objective function J(z ) is given by 



"W-2 E,(zmz) dT( *> 



(3.29) 



dkj 4 v~/' * v~y dk _ 

Thus in order to find the gradient with respect to the adaptive pa- 
rameters k : , we must take the partial derivative of T(z ) with respect to each kj 



ST(z) 

dkj 



i=i 

i*j 



dkj 



(3.30) 
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Kwan and Martin Structure 




Fig 3.12 



From equation (3.24) we have 



dHjAz) = d[l - H 3 bp (z) 1 
Ok, dk : 

d H^z) 

dk 3 



(3.31) 

(3.32) 



Using the rule for differentiation as given by the equations (2.45) and (2.47) and from 
equation (3.23) we have 



8uk*) _ 

D(z) 



(3.33) 

(3.34) 



where D(z ) is the denominator of Hl p . H{, is given by the equation (3.15). 

Substituting equation (3.34) into equation (3.30) we obtain the over 
all sensitivity for the jth ki parameter as 



^ = i\ir x (z)H’ tr (z)H’,(z) 



1=1 



fffli-M 

t=l 



n 

»=j-i 



Hl r (z)H{,(z) 



By recognizing that 



A T 



n hu*)k 

«=j- i 



Y’{z) = 

and by substituting equation (3.37) in the equation (3.36) we get 

dT(z) 



(3.35) 

(3.36) 

(3.37) 



dkj 



n^(-) 

i=i 



Hi,(z) 



(3.38) 



Figure 3.12 shows the Kwan and Martin realization of the complete 
adaptive system. The difficulty in the Kwan and Martin approach is in the generating 
of the product of notch filters without the notch filter j, as required in equation (3.36). 
To generate this product for each section, w'ould require JV — 1 biquads per section 
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resulting in a total of N(N — 1) biquads just to generate the product. Kwan and 
Martin are able to reduce this by using the output of the bandpass filter as the input 
to their cascade via equation (3.37). Since this output already has ( j — 1) of the 
required H' N (z) factors in it, only (N — j) additional biquads are needed for a total of 
0.5(A r2 — N). Adding this to the N biquads required to realize the cascade of notch 
filters and the N biquads required to realize the H : s (z) factors, yields a total number 
of biquads given as 0.5(A r + 3)A r . 

2. The New Structure 

Figure 3.13 shows the improved adaptive notch filter structure [Ref. 25], [Ref. 
26]. The key to the improvement is the recognition that the output Ei(z) = T(z)X(z) 
for the cascade of the notch filters can be written both as a product of the individual 
notch filter section transfer functions H' N (z ) times the input and in terms of the input 
X(z) minus the outputs Y'(z ) of the bandpass filters 



Er(z) = T(z)X(z) 

N 



n^ivM 



i=l 



X(z) 



= A'(z)-(£r( 2 )) 



(3.39) 

(3.40) 

(3.41) 



To get the product of H' N (z) without the term i = j as required in the 
equation (3.36), we may use equation (3.41) to simply add back the term Y*(z) 



rW(-) 


A'(z) 


= X(z)~ 


N 


is=l 






»=1 



(3.42) 

= E x (z) + Y 3 {z) (3.43) 

Figure 3.13 makes use of this fact to generate the gradient needed for the 
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adaptive process. From the Figure 3.13, we can see that the total number of biquads 
required is N for the cascade of notch filters plus 2N for the H bp (z)H' t (z ) required for 
adaptation, minus 1 at the last stage, since the last stage does not need the extra 
H' bp (z)). Thus we have the number of biquads required in the new structure is (3N-1). 
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IV. MODELLING and SIMULATION 



A. MODELLING OF SIGNALS 

In order to test various algorithms and to evaluate their probable performance 
in the real environment, it was necessary to develop a meaningful simulation of the 
real situation. To achieve this, the following four testing categories were developed: 

(i) Sinusoidal signals with white gaussian noise 

(ii) Narrow Band Noise with white gaussian noise 

(iii) Bi-Phase Shift Keying (BPSK) sequence 

(iv) Frequency Shift Keying (FSK) sequence 

1. Sinusoidal Signals 

In order to generate sinusoidal signals with minimum computational burden 
placed at different frequencies a second order AR process 

x\ = 2cos(9i)x' k _ 1 - x\_ 2 (4.1) 

was used. Initial conditions are very important and they are chosen such that x_i = 0 
and x _2 = —sin(O t ) giving a unit amplitude sinusoidal signal. The # t value is between 
0 to 180 and n is the number of frequencies desired. The required signal y * needed 
to input into the adaptive algorithm is given as 

n 

y* = H** + 7* (4.2) 

1=1 

where 7 * is a white gaussian noise N(0,a 2 ). 

2. Narrow Band Noise 

The narrow band noise signal is generated using the difference equation 
4 = 2rcos(0,)4_, - r 2 x \_ 2 -f u\ (4.3) 
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where 0, decides the placement of the noise in the spectrum and r controls the band- 
width of the noise. The u\ is simply a uniformly distributed noise taken at different 
instants. The desired signal yk is obtained via 

Vk = 'Y,z' k + lk (4.4) 

«=i 

3. Bi-Phase Shift Keying Sequence 

The generation of BPSK signal has three distinct three parts: 

(i) Generation of Random Binary Sequence(RBS) is achieved by passing 
uniformly distributed noise through a hardlimiter (An important note is that the 
interval between the two consecutive bits of RBS is j^); 

(ii) Another sequence of binary numbers is a spreading code or sequence. 
The specific sequence used in a given communications system is normally not available 
to anyone but to the designated receiver. (In this particular simulation we have 
generated the spreading sequence by passing uniformly distributed noise through a 
hardlimiter. The bit interval is 7-); 

(iii) Phase encoding ( i.e .) mapping the given binary signal w-hich is the 
exclusive or of (i) and (ii) as 0 or tt at appropriate sample time. 

The output of the first hardlimiter is stored in an array x. Output of the 
second hardlimiter is stored in the array y This information is retrieved by a subtle 
use of the array index given as i = k /(, where i is the index, k is the discrete sample 
number and /<, is the bit rate of the intelligence. Similarly another index j is generated 
using j = kf c where f c is the chip frequency. Generation of the indices is the key 
thing in this simulation. The desired signal now is given by the equations 



wl = Acos(2n f 0 k -f <f>k) 


(4.5) 


<t>k = [*(*’) © y(j)]7r 


(4.6) 
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a k = ^2w p k (4.7) 

p=l 

where i and j are indices of the arrays as defined earlier. It is an important point 
to note that the characteristic of the signal w p k are defined by the parameter p = 
{ foi /c? fb} • 

This signal is not really a simple BPSK signal but it has an additional fea- 
ture of spreading the spectrum by controlling the chip-frequency and carrier frequency 
and information rate. A block diagram of the scheme is given in Fig 4.1. 

The desired signal y k is given by: 

y k = a k + fi k (4.8) 

where a k is the set of narrowband BPSK signals placed at different places in the 
frequency spectrum and 0 k is the broad band BPSK signal generated for a specific p 
value. 

4. Frequency Shift Keying Sequence 

Generating this sequence needs a random binary intelligence signal. This 
was once again is achieved by passing a uniformly distributed noise through a hardlim- 
iter. The output of this hard limiter stored in an array x. An index i is chosen such 
that i = k /j, where fo is the baud-rate of the information and k is discrete sample 
number. Now the desired signal is generated via 

s k = 2cos(0*)s*_i - s k . 2 (4.9) 

0 k = 0 + 6x(i) (4.10) 

Vk = s k + 7 * (4.11) 

where 6 is the carrier frequency and <5 is the depth of the frequency modulation. 
Initial conditions are very important and they are chosen such that s_i = 0 and 
2 = —sin{0) giving an unit amplitude sinusoidal signal. 
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B. SIMULATION 



The adaptive digital filter algorithm as described in the earlier chapters is sim- 
ulated and tested using synthetic data. Simulation was carried out on a VAX 11/785 
computer using Fortran 77. Listing of the program used is enclosed in the appendix. 
The adaptive filter parameters of interest are 

(a) Sharpness of the notch filter defined by pole position (r 2 = 1 — ac 2 ) 

(b) Step size in the parameter update procedure (fi) 

(c) Time constant of the fading filter (A) ( refer to equation 3.19) 

(d) Model order (N= # of 2nd order filters) 

(e) Order of the incoming signal or number of interferers (m) 



1. Kwan Martin Algorithm 

In simulating this algorithm, the most important thing is the implemen- 
tation of the notch filter. Fig 4.2 gives the structure of the algorithm. Let x(n) be 
the input to the Adaptive Filter and e(n) the desired output. This desired output is 
obtained by passing the input x(n) via a cascade of N notch filters as 

e(n) = {//"(r)//#- 1 • • • H' N (z)} x(n) ' (4.12) 

Output at the intermediate jth section of the band-pass filter is designated as y*(n) 
as shown in Fig 4.2 Then the sensitivity of the parameter, is obtained by passing 
this y J '(n) through another cascade of (j — 1) notch filters given as 

s‘(n) = y>(n) (4.13) 

The filter transfer functions H J N (z ) and H> 3 were defined in the earlier chapter by 
equations (2.21) and (2.15). 
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2. The New Algorithm 

The primary difference between the New Algorithm and the Kwan and 
Martin algorithm is in the method of obtaining the sensitivity of the j th coefficient. 
This difference can be seen in the Fig 4.3. The output t/ J (n) at the jth notch Fig 4.3. 
filter is added to the over-all output e{n) and this summed output is passed through 
a cascade of only two filters as shown in Fig 4.3 and can be given as 

s’(n) = [h^hQ {e(n) + y’(n)) (4.14) 



3. Forgetting Filters 

For parameter incrementation via equation (3.17) we need to obtain ||s 2 (n)||. 
This could be achieved by 

a) Zero order forgetting using equation (3.22) 

b) 1st order forgetting using equation (3.19) 

c) 2nd order forgetting using equation (3.21) 

4. Stability 

The parameter incrementation given by the equation (3.17) can be recast 
by using either of the forgetting filters given above as 



k\ = k\ — //[*/] 1 s'(n)e(n) 



(4.15) 



It is very important to note that the above equation (4.15) has a close resemblance 
with the ML Estimator. In the MLE case the stability while incrementing the param- 
eter is an important factor. A similar problem exists in the current incrementation 
procedure but the problem is solved by a suitable choice of parameters and modify- 
ing the incrementation procedure by adding an additional factor to equation (4.15) 
as follows: 



k\ = k\ — fie(n)s'(n) 



1 + 



Vma 



Vmin + V 



(4.16) 
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This modification also protects against possible overflow or underflow while computing 
[v‘] . In addition simulations indicate that in the case of zero-order forgetting, this 

factor is required in order to obtain convergence. When u mtn and u max are zero and 
infinity respectively, equation (4.16) reduces to equation (4.15). 

In this incrementation procedure two forms of the denominator polynomial 
of the Hb p (z) were considered: 



D'{z) = 1 - (2 -k 2 - ki)z~ l + (1 - k 2 )z 



.-2 



(4.17) 



and 



D'(z) = 1 -2 K\rz-' +r 2 z~ 2 



(4.18) 



where k\ = cos(O'). 

While using the incrementation procedure, under transient conditions poles 
of the D'(z) cross the unit circle causing instability problems. This condition was 
averted by checking the pole position after incrementation and if unstable then in- 
crementation is modified. Checking this condition for D'{z ) given by equation (4.17) 
calls for solving a quadratic equation at every parameter increment and examining 
the pole position. However this check is very simple for D'(z) given in equation 
(4.18), since we need only to check by maintaining |k\| < 1. Furthermore the 
value of r must be positive and less than unity for stability. The choice of the r is 
very important for proper fast convergence. Fig 4.4. shows the frequency response of 
band-pass filter for different values of r The 3dB band-width of the band pass filter 
is related to r and is given [Ref. 29] by 



Band — width = cos 




(4.19) 



Note that under limiting conditions, the band-width is essentially zero. In this simu- 
lation it was assumed a value of r as 0.95. 
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A large number of simulations [Ref. 29] were carried out using sinusoidal 
inputs with and without noise for different m and N values. In these simulations 
the denominator polynomial used was given by equation (4.18). The results were 
tabulated in Table 3.1. From Table 3.1 it is seen that the choice of the step size is 
an important factor and step-size must be optimized for a specific number of sections 
N. The values of u mt -„ and v max did not pose any problem while incrementing using 
1st order or 2nd order forgetting filters. It seems that a 1st order forgetting filter is 
more effective than either zero-order or a 2nd order filter. 

After fixing the values of fi and A the algorithm was tested using Narrow 
Band Noise signal y* for its performance. This signal was used only to tune the value 
of r (i.e. sharpness of the notch filter). The following simulation results are based on 
parameter adaptation for the D'(z ) given by equation (4.17). 

5. Response for Sinusoidal signal 

A sample simulation output due to sinusoidal signal was shown in Fig [4.5 
- 4.7]. The input signal is composed of 3 sinusoids with normalized frequencies 
JUj /,, and Hc»/j aQ d WGN with <7 = 1. The spectrum of this signal is shown in Fig 
4.5. After passing this signal through the adaptive filter we could see that interferers 
were removed and only the noise was left behind. This is clearly shown in Fig 4.6. 
The adaptation process is shown in Fig 4.7. 

6. Response for Narrow Band Noise 

In this testing category sinusoids were replaced by narrow-band signals. 
These signals were generated via equation (4.4). Superimposed on this signal, WGN 
with <7 = 1 was added. A typical spectrum of this signal is shown in Fig 4.8. This 
signal was passed through the adaptive filter and these narrow band interferers are 
notched out and only WGN is left out as showm in the Fig 4.9. Fig 4.10 shows the 
corresponding parameter convergence. 
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interferers 0.95 

Variance of the noise 1. 

Adaptive filter step . . . .0.01 
Adaptive filter r 0.9 



7. Tracking FSK Signal 

Transient behavior of the adaptive filter was studied by applying an FSK 
signal (generated via equation 4.11). In the case of an FSK signal, signal spectrum 
constantly changes. This property is useful for testing the tracking ability of the 
adaptive filter. In fact a WGN with a = 1 was also added to the signal. An adaptive 
filter then was used to demodulate the signal, by tracking the spectrum. This is 
shown in Fig 4.13. 

8. Suppressing BPSK Interference 

Having seen the performance of the adaptive filter under various conditions, 
it is only needed to test under an additional constraint of broad band noise. In this 
case, the broad-band signal was swamped by narrow-band interferers and WGN. A 
typical broad-band signal is generated by using a BPSK signal generator via equation 
(4.8). In this we have used carrier at 0.25 cycle/sample and the chip frequency at 0.125 
cycles/sample. The narrow-band interference signals are also generated by using the 
same BPSK signal generator, but with different parameters. Interferers are chosen 
such that the chip frequency is fixed at ^ f s cycles/sample, while carriers were chosen 
at 360 ^*’ §6of‘i an ^ 36 o-^ s cycles/sample. To these interferers and broad-band signal, 
a WGN with <7 = 1 was added. The composite signal spectrum is shown in Fig 4.14, 
indicating the three interferers, broad-band signal, and WGN. This signal was passed 
through the adaptive filter. Output of the filter is shown in Fig 4.15. In this figure 
it is clear that the interferers were removed by the filter and only the desired signal 
was left behind. 

9. Model order mismatch 

The model order of the algorithm was fixed at 3 while BPSK signals were 
generated with 2 interferers i.e. m = 2 and N = 3. Fig 3.16 shows the parameter con- 
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vergence under these conditions and Fig 3.19 shows for m = 1 and TV = 3 conditions. 
We could notice that free extra parameter(s) wandering due to the mismatch. This 
has the effect of notching out the desired signal. This can be solved by deliberately 
introducing a known BPSK signal at fixed place. 
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Adaptive filter 



Summary of Tests on the New Algorithm 



TABLE 4.1 



Case 


# of 
sections 


Sines 


^mtn 


^max 


Noise 


zero order 


1st order 


2nd order 




iter 




iter 


t* 


iter 


la 


i 


i 


0.05 


20 


no 


0.05 


67 


0.05 


20 


0.05 


20 


lan 


i 


i 


0.05 


20 


yes 


0.05 


jump 


0.05 


20 


0.05 


144 


lb 


i 


3 


0.05 


20 


no 


0.05 


jump 


0.05 


30 


0.05 


25 


lbn 


i 


3 


0.05 


20 


yes 


0.05 


jump 


0.05 


50 


0.05 


100 


2a 


5 


3 


0.10 


10 


no 


0.017 


340 


0.063 


200 


0.07 


400 


2an 


5 


3 


0.10 


10 


yes 


0.017 


1200 


0.063 


400 


0.07 


450 


2b 


5 


5 


0.10 


10 


no 


0.017 


1410 


0.063 


150 


0.07 


500 


2bn 


5 


5 


0.10 


10 


yes 


0.017 


3500 


0.063 


175 


0.07 


250 


2c 


5 


7 


0.10 


10 


no 


0.017 


1667 


0.063 


450 


0.07 


jump 


2cn 


5 


7 


0.10 


10 


yes 


0.017 


700 


0.063 


167 


0.07 


jump 


3a 


10 


7 


0.02 


50 


no 


0.02 


1667 


0.075 


500 


0.07 


1500 


3 an 


10 


7 


0.02 


50 


yes 


0.02 


5000 


0.075 


750 


0.07 


1000 


3b 


10 


10 


0.02 


50 


no 


0.02 


7500 


0.075 


2000 


0.07 


5250 


3bn 


10 


10 


0.02 


50 


yes 


0.02 


1000 


0.075 


2000 


0.07 


2500 



• In all cases with noise, parameters converge to exact value when number of 



notches is greater or equal to number of sinusoids. Iterations listed in table 
represent convergence to three decimal places. 

• In all cases with noise, parameter oscillate around the exact value. The number 
of iterations listed in the table is the number of iterations required to establish 
this oscillatory pattern. 

• In some cases with more interfering sine waves than notches will jump at random 
between sine waves. This pattern is indicated in the table by the word jump. 

• Values of u m ,- n and v max in the table refer only to zero-order forgetting which 
generally will not converge without limits on v. All other types' of forgetting 
were run without limits on v 
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V. HARDWARE IMPLEMENTATION 



Although hardware design is beyond the scope of this thesis, in this chapter we 
shall briefly look at possible hardware configurations with a view toward feasibility 
of the new algorithm. 

A. HARDWARE FEASIBILITY 

In recent years many dedicated digital signal processing chips have become 
commercially available. Table 5.1 gives characteristics of some typical chips that are 
currently available. 



TABLE 5.1 



Feature 


Commercial Make 


MIPS R3010 


Weitek 3364 


TI 8847 


AT&T DSP32c' 


Cycle Time(ns) 


40 


50 


30 


20 


Cycles/add 


2 


2 


2 


2 


Cycles/mult 


5 


2 


3 


2 


Cycles/divide 


19 


17 


11 


3 


Cycles/sq root 


- 


30 


14 


? 



Data on MIPS R3010, Weitek 3364, and TI 8847 was obtained from the computer 



architecture book by Hennesey and Patetrson [Ref. 30] while the data for AT&:T 
DSP32c is obtained from DSP32c data manual [Ref. 31]. Fig 5.1. gives a schematic 
of the 2nd order HR filter hardware scheme. It has a coefficient memory which can 
be set by an external device. Close examination reveals that Multiplications A can 
be performed simultaneously while Multiplication B and other additions are to be in 
sequence. Fig 5.2. is the basic building block for the adaptive filter and is identical for 
all the sections. This is offers an advantage in hardware complexity over the Kwan <*c 
Martin approach. Filters 1 and 2 have transfer functions of Hb p (z) while filter 3 has a 



Ox ka{x 4c*k T> do ‘a • 
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transfer function H 3 ,(z). Throughput rate is primarily determined by the computing 
time of the 2nd order HR filter. The overall block diagram of the adaptive filter for a 
3-Notch cancellation was given in Fig 5.3. This architecture remains same either for 
a Floating point or Fixed Point arithmetic. 

1. Time Budget 

In the architecture of Fig 5.1, 5.2, and 5.3 the most vital elements are the 
IIR filter , Controller and the Forgetting filter. The HR filter corresponds to Hb p {z). 
The controller has to update the parameter via equation (3.17) which calls for 2 
multiplications, 1 addition, and 1 division. Similarly the Forgetting Filter has to 
implement equation (3.19) calling for 2 multiplications and 1 addition. Results are 
tabulated as given below 



TABLE 5.2 



Element 


# of Effective 
Multiplications 


# of Effective 
Additions 


# of 
Divisions 


HR filter 


2 


3 


nil 


Controller 


2 


1 


i 


Fogetting 

Filter 


2 


1 


nil 



The maximum computing time is at the IIR filter. Using the AT&T DSP32c proces- 
sor, it takes 5x2 = 10 cylces (ie 200ns) for implementation the desired IIR filter. This 
corresponds to a Throughput rate of 5 Mhz. This is the best that could be achieved 
with the existing commercial floating point processors. However processors like the 
HDSP66110 of Honey well make claim an ALU speed of 10ns with block floating point 
operations can conveniently implement at 10 Mhz throughput rate, without pipelin- 
ing the IIR filter. With pipelining the throughput rate would increase dramatically, 
but will never exceed the maximum throughput rate of 100 Mhz. 
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2. VLSI Approach 

Taking advantage of the fact that the entire hardware can be generated by 
a simple repetitive use of the building block, strong consideration should be given 
to designing a VLSI chip for Fig 5.2 rather than implementing the hardware via the 
existing commercial chips. The main rationale behind this idea is that reliability 
of the system is very important in EW equipment. Also VLSI has the potential of 
ultimate-low cost. 
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VI. CONCLUSIONS 



In this thesis the problem of suppressing narrow-band interference was ad- 
dressed. The problem specific to the Electronic warfare scenario was kept in mind 
while solving the problem. The new algorithm [Ref. 26] was derived and simulated. 
This new algorithm was tested against various signals and signal conditions. The 
results are highly encouraging. 

Some aspects of pipelining were discussed. Reduced hardware complexity of the 
New Algorithm is an important advantage over the earlier algorithm by Kwan and 
Martin [Ref. 8]. A possible hardware scheme for implementing the new algorithm 
was discussed. It was observed that with the existing commercially available floating 
point processors a throughput rate of 5 Mhz is achievable while using processors like 
HDSP66110 with an ALU speed of 10ns it is possible to achieve a throughput rate 
of 10 Mhz. 

A. FUTURE WORK 

Future directions of work are 

i) VLSI design of the system with an architecture 
as indicated earlier or a similar architecture, 

ii) Pipelining of 2nd Order HR filter suited for 
this application 

iii) Design using recursive DFT 

The above areas of w’ork are the logical extensions of this work. However a radically 
different approach for the same problem using adaptive filtering in frequency domain 
[Ref. 10] should also be considered. 
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APPENDIX A 



dimension faray(5000) 
dimension uaray(5000) ,yaray(5000) 
dimension fhat (5000) ,ph (5000) 
dimension a(90),b(90) 

open (unit=9, f ile= ' spk. dat ' , status='new' ) 
c 

c This programme generates input output 
c sequence by exciting a linear system 
c defined by the numerator polynomial 
c B(z) and denominator polynomial A(z) 
c this data is stored in uarray and yarray 
c Subsequently an ARMA (pole, zero) is used 
c to fit this data, 

c 

n=1024 

ix=l 

yk=rand ( ix) 
ix=0 
kk=8 
c 

c generate sequence uk and yk 

c 

ip=3 

iz=2 

a (1) =1. 7 
a ( 2 ) =-l . 53 
a ( 3 ) =0 . 648 
b (1) =1 . 0 
b (2 ) =0 . 6 
do 10 k=l,n 
call rbn(uk,ix,k) 
call system(a,b,uk,yk, ip, iz, k) 
uaray (k) =uk 
yaray (k)=yk 
10 continue 
c 

c save the sequence in the arrays uaray and yaray 
c 

call spktrm (yaray , faray) 
c 

c spectrum of the output sequence yaray is computed 
c using FFT programme and stored in faray for 

c comparision 

c 

ip=4 

iz=l 

c 

c Fit an ARMA ( ip , iz ) model for the given data 

c this is an ordinary least squares algorithm 
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c 

print * , ' input-ip-iz — >' 
read *,ip,iz 
do 20 k=l , n 
uk=uaray (k) 
yk=yaray (k) 

call anna (a ,b, uk,yk, ip, iz , kk, k) 
c print * , ' num-> • , (b ( i) , i=l , iz) 
c print * , 'den-> ' , (a (i) , i=l , ip) 

20 continue 

print * , ' num-> • , (b ( i) , i=l , iz) 
print * , 'den-> ' , (a (i) , i=l , ip) 
c 

c Sampling of the Z-transform. given 

c coefficients a .... and b 

c corresponding H(z) = B(z)/A(z) 

c is evaluated on the unit circle 

c for obtaing the parametric spectrum 

c 

call zsmple (a,b, fhat,ph, ip, iz) 
c 
c 

do 30 i=l,n/2 

c print *,i,' — > ' , fhat (i) , faray ( i) 

write (9,*) i, fhat(i) , faray (i) 

30 continue 

stop 
end 
c 
c 

subroutine zsmple (a , b, r , ph, ip, iz) 

dimension a(90),b(90) 

dimension r (2048) ,ph(2048) 

complex z , ui , fn, fd , omega , delw, spec 

pi=atan (1.0) 

pi=4 . 0*pi 

ui=(0.0, 1.0) 

delw= (pi/1024 .0) *ui*2 . 0 

omega= (0 . 0 , 0 . 0) 

do 100 k=l ,1024 

z=exp (-omega) 

fd=(0. 0,0.0) 

do 10 i=l,ip 

10 fd=fd+a(i)*(z**(-i) ) 

fd=l . 0-fd 
fn= (0 .0,0.0) 
do 20 i=l,iz 

20 fn=fn+b(i) * (z** (-i) ) 
spec=fn/fd 
r (k) =abs (spec) 
specr=real (spec) 
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speci=aimag (spec) 
ratio=speci/specr 
ph (k) =atand (ratio) 
omega=omega+delw 
c print — z-om-> z , omega 

100 continue 
return 
end 
c 
c 

function atand(x) 
pi=atan (1.0) 
pi=4 . 0*pi 
xrad=atan (x) 
atand=xrad* ( 180 . 0/pi) 
return 
end 
c 
c 

subroutine rbn(b,ix,k) 
z=rand ( ix) -0 . 5 
if(z.ge.O.O) b=1.0 
if(z.lt.O.O) b=-l . 0 
return 
end 

c 

c 

subroutine system(a,b,uk,yk, ip, iz,k) 
dimension 

& zkar (90) , zkma(90) , 

& a(90) ,b(90) 

if (k.ne.l) go to 250 
do 230 i=l,ip 
230 zkar(i)=0.0 

do 240 i=l,iz 
240 zkma(i)=0.0 

250 continue 

do 30 i=ip,2,-l 

30 zkar ( i) =zkar ( i-1) 
zkar (l)=yk 

do 31 i=iz,2,-l 

31 zkma (i) =zkma (i-l) 
zkma (1) =uk 

c 

c 

yk=0 . 0 

do 200 i=l,ip 
200 yk=yk+a ( i) *zkar (i) 
do 210 i=l,iz 
210 yk=yk+zkma (i) *b (i) 
return 
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end 

c 

c 

subroutine delay (uk, yk, idelay, k) 
dimension d(500) 
if (k.ne.l) goto 10 
do 20 i=l , 500 
20 d ( i) =0 . 0 
10 continue 

do 30 i=500 , 2 , — 1 
30 d ( i) =d ( i-1) 
d ( 1) =uk 
yk=d (idelay+1) 
return 
end 
c 
c 

subroutine spktrm(taray , faray) 
complex u(5000) 

dimension taray(5000) , faray (5000) 

m=10 

n=1024 

rn=n 

do 20 i=l,n 
u ( i) =taray ( i) 

20 continue 

call pfft(u,m,n) 
do 30 i=l,n/2 
faray (i) =real (u ( i) ) 

30 continue 

return 
end 
c 
c 

c ********************************** 

c * 

c Subroutine for computing DFT of 

c an array 'a' is complex and a pair 

c of numbers are to be specified 

c for each point * 

c m is the 2 power index * 

c say m=10 then number of points * 

c are 1024 * 

c after computation fft is kept * 

c in the same complex array 'a' * 

c * 

0 ********************************** 

c 

subroutine pfft(a,m,n) 
complex a (5000) ,u,w, t 
complex ui,ur 
pi=3. 141592653589793 
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n=2**m 

nv2=n/2 

nml=n-l 

j=l 

do 10 i=l,nml 
if(i.ge.j) goto 20 
t=a ( j ) 
a( j)=a(i) 
a ( i) =t 
20 k=nv2 

30 if(k.ge.j) goto 10 
j=j-k 
k=k/2 
goto 30 
10 j=j+k 

do 40 1=1, m 
le=2**l 
lel=le/2 
u=(l. 0,0.0) 

w=cmplx (cos (pi/lel) , -sin (pi/lel) ) 
do 40 j=l,lel 
do 50 i=j,n,le 
ip=i+lel 
t=a ( ip) *u 
a ( ip) =a (i) -t 
50 a(i)=a(i)+t 

40 u=u*w 
rn=n/2 

ui= (0 . 0,1.0) 
ur= ( 1 . 0 , 0 . 0) 
do 60 i=l,n/2 
a(i)=a(i)/rn 
a ( i+512 ) =a ( i+512 ) /rn 
areal=real (a (i) ) 
a imgn=a imag ( a ( i ) ) 
amagn=abs (a ( i) ) 
aphase=atand ( aimgn/areal ) 
a ( i ) =ur*amagn+ui*aphase 
areal=real (a (i+512 ) ) 
aimgn=aimag(a(i+512) ) 
amagn=abs ( a ( i+512 ) ) 
aphase=atand (aimgn/areal ) 
a ( i+512 ) =ur*amagn+ui*aphase 
60 continue 
return 
end 
c 
c 

c $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 

c 

c 

subroutine anna (a , b, uk, yk, ip, iz , kk, k) 
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dimension 

& zkar (90) , zkma (90) , 

& zk(90) ,a(90) ,b(90) ,c(90) 

if (k.ne.l) go to 250 
n=ip+iz 
ipp=ip+l 
do 230 i=l,ip 
230 a (i) =0 . 1 

do 240 i=l,iz 
240 b ( i) =0 . 1 
250 continue 

do 30 i=ip,2,-l 

30 zkar(i)=zkar (i-1) 
zkar (l)=ykold 
ykold=yk 

do 31 i=iz,2,-l 

31 zkma (i) =zkma (i-1) 
zkma (1) =uk 

do 32 i=l,ip 
c(i)=a(i) 

32 zk(i)=zkar (i) 
do 33 i=l,iz 
c ( i+ip) =b (i) 

33 zk(i+ip) =zkma ( i) 

call syseq (c , yk, zk, n , kk, k) 
do 34 i=l,ip 

34 a(i)=c(i) 

do 35 i=l,iz 

35 b(i)=c(i+ip) 

return 

end 

c 

c 

c 

subroutine syseq (a , yk, zk, n, kk, k) 
c 

c *************************************** 

c solves system of equations with 

c # of equations more than unknowns 

c using linear reggression. 

c t 

c equation is yk=zk * a 

c where 'a' is n vector to be estimated, 

c 

c **************************************** 

c 

dimension 

& del (90) , phat (90,90) , zk (90) , y (90) , a (90) 
if (k.ne.l) go to 250 
alpha=l . 0 
250 continue 

call inv (phat , zk, alpha , n, kk, k) 
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ykhat=0 . 0 
do 40 i=l,n 

40 ykhat=ykhat+zk ( i) *a ( i) 

er=yk-ykhat 
do 60 i=l,n 
60 y (i) =er*zk(i) 

do 80 i=l,n 
del ( i) =0 . 0 
do 80 j=l,n 

80 del (i)=phat(i, j ) *y ( j )+del (i) 

if (mod(k, kk) .ne. 0) goto 70 
c print * , • > ' , k 

c print 100, ( (phat(i, j) , i=l,n) , j=l,n) 

c print 130, (del (i) , i=l, n) 

c print 160, er 

c print 170, (zk( j ) , j=l , n) 

70 continue 

do 50 i=l,n 
50 a (i) =a ( i) +del (i) 

100 format (2x, 'phat ', 2x,/ , 4el6. 6) 

120 format(2x, ' zk' , 2x,/,4el6.6) 

130 format (2x, 'del ' , 2x,/, 4el6. 6) 

160 format(2x, ' ekl ' ,2x,el6.6,/) 

170 format (2x, ' zk ', 2x,/ , 4el6 . 6) 
return 

end 

c 

c 

subroutine inv (phat , zk, alpha, n, kk, k) 

dimension phat(90,90) ,zk(90) ,q(90) ,r(90) ,delp(90,90) 
if(k.ne.l) goto 5 
do 6 i=l,n 
do 6 j=l,n 

if(i.eq.j) phat (i , j ) =1 . 0 
if(i.ne.j) phat ( i , j ) =0 . 0 
6 continue 

5 continue 

do 90 i=l,n 

do 90 j=l,n 

phat ( i , j ) =phat ( i , j ) /alpha 
90 continue 

do 10 i=l,n 

q ( i) =0 . 0 
do 10 j=l,n 

10 q(i)=phat(i, j) *zk(j)+q(i) 
sum=l . 0 

do 20 i=l,n 

20 sum=sum+zk(i) *q(i) 
do 30 j=l,n 
r(j)=0.0 
do 30 i=l,n 

30 r ( j ) =phat ( i, j ) *zk (i) +r ( j ) 
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do 40 j=l,n 
do 40 i=l,n 
delp(i, j)=q(i) *r( j) 
do 50 j=l,n 
do 50 i=l,n 
50 phat (i, j ) =phat (i , j ) - (delp (i, j )/sum) 

return 
end 
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c 

c 

c 



APPENDIX B 



simulation of michal paper IEEE 

dimension akl (10) , ak2 (10) , s (10) , xout (10) , theta (10) 
dimension aks(10) ,wp(10) ,ak3 (10) 
dimension array (4,4000) 

open( unit=9, file=' martin.dat' , status= 1 new ' ) 
open( unit=8, file='mart.dat ' , status= 1 new' ) 

print * , ' -input-SNR-in-dB >' 

read *,snr 
c if (snr.ne.0) goto 200 

c var=0.5 
c goto 210 

var=l .0/(2.0*(10.0** (snr/10. 0) ) ) 

210 continue 

print — variance — >',var 

n=l 

print *,'-input-nh — >' 
read *,n 

print * , ' -input-#-of-waves — >' 
read *,nd 

print *,' -input-angles — >' 
read *, (theta (i) , i=l,n) 
r=0 . 9 

pi=4 . 0*atan (1.0) 
rad=pi/180 . 0 
amu=0 . 0 

c do 50 kl=l , 50 
amu=0 . 015 

print -input-step-size — 0.015 — >',amu 

read * , amu 

print *, ' — angles — >' , (theta (i) , i=l,n) 

print *, ' — input — 1-for-step-change — else-0 — >' 

read *,ic 

if (ic.eq.l) print *,' — input-step-in-degrees — >' 

if (ic.eq.l) read *,step 

do 40 i=l,n 

ak2 (i)=(l-r*r) 

ux=theta ( i ) *rad 

value=cos (ux) 

c akl ( i) =dl*sqrt ( 1 . 0+r*r-2 . 0*r*value) 

akl ( i) =sqrt ( 1 . 0+r*r) 
aks (i) =sqrt ( 1 . 0+r*r-2 . 0*r*value) 

40 continue 
c amu=amu+0 . 0001 

do 10 k=l ,6500 

call system(yk,wk, theta, step, var, ic,nd,k) 
call f il ters ( akl , ak2 , aks , s , en , xout , yk , n , k) 
call increment ( akl , ak2 , aks , s , en , amu , avg , n , k) 
do 30 i=l,n 
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xx=2 . 0*sqrt (1-0. 5*ak2 (i) ) 
wp ( i) =2 . 0*asin (akl (i) /xx) 
c wp(i)=wp(i)/rad 
wp(i)=cos(wp(i) ) 

30 continue 

do 60 i=l,n 

ak3 (i) =abs (akl ( i) ) 

60 continue 

kk=mod (k, 100) 
ll=mod (k, 10) 

if(kk.eq.O) print *, ' — akl — >' , (ak3 ( j) , j=l,n) 
if(ll.eq.O) write(9,*) k, (ak3 ( j ) , j=l , n) 
if(kk.eq.O) print ' — aks — > ' , (aks ( j ) , j=l, n) 
c if(kk.eq.O) print *,k, ' -avg->' , avg 
c print * , (wp(i) , i=l, n) 

c write(9,*) k, (wp ( i) , i=l , n) 

if (k.gt . 1500. and.k. It . 2500) write(8,*) k,yk,en 
if (k.gt . 5000. and. k. It . 6000) write(8,*) k,yk,en 
10 continue 

c print *, kl,amu,avg 

c write (9,*) kl,amu,avg 

50 continue 

close (9) 
stop 
end 
c 
c 

subroutine increment ( akl , ak2 , aks , s , en , amu , avg , n , k) 
dimension akl (10) ,ak2 (10) , aks (10) , s(10) , ss(10) 
dimension dec (10) , flag (10) , fading (10) 
if (k.ne.l) goto 30 
c amu=0.001 

avg=0 . 0 
fading ( 1) =0 . 9 
fading (2) =0.9 
fading ( 3 ) =0 . 9 
fading(4)=0.9 
30 continue 

do 10 i=l,n 

ss (i)=fading(i) *ss(i)+s(i)*s(i)*(l. 0-fading (i) ) 

10 continue 

c do 31 i=l,n 

c3 1 print * , ' -en-s-ss-> ' , en, s (i) , ss (i) 
c sum=0 . 0 

do 20 i=l,n 

dec ( i) =amu*en*s (i)/(ss(i)+0.0001) 
akl ( i ) =akl ( i ) -dec ( i ) 

c sum=sum+ (abs (aks ( i) -akl ( i) ) )/ (aks ( i) ) 

call stability (akl, ak2, flag, n) 
akl (i) =akl ( i) +flag ( i) *dec (i) 

20 continue 

realk=k 
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avg=avg* ( (realk-1 . 0) /realk) +sum/realk 
return 
end 
c 
c 
c 

subroutine filters ( akl , ak2 , aks , s , en , xout , xin , n , ki ) 
dimension w(10,3) , a(10,2) , b(10,2) ,e(ll) , xbp(10) 
dimension g(10, 3) ,gbp(10) ,u(10) 
dimension q(10,3) ,s(10) ,xout(10) ,c(10) 
dimension akl (10) , ak2 (10) , aks (10) , as (10) 
c 

if (ki.ne.l) goto 200 

c open( unit=9 , file= 'martin.dat 1 , status= 'new' ) 

200 continue 
i=n 

do 40 i=l,n 

a (i , 1) =- (akl ( i) *akl (i) +ak2 ( i) -2 . 0) 
as ( i) =- (aks ( i) *aks ( i) +ak2 ( i) -2 . 0) 
a (i , 2 ) =- ( 1 . 0-ak2 ( i) ) 
b(i , 1) =-0 . 5*ak2 (i) 
b ( i , 2 ) =0 . 5*ak2 ( i) 
c(i) =-akl (i) *ak2 (i) 

40 continue 
c do 41 i=l,n 

c print * , ' -kl-k2-> ' , akl (i) ,ak2 (i) 

c41 print * , '-a-b-c-> ' , a (i , 1) , a (i, 2) ,b (i, 1) ,b(i , 2) ,c (i) 

c 

c 

e(l)=xin 

i=n 

do 10 j=l,n 
w(i/3)=w(i,2) 
w ( i t 2 ) =w ( i f i ) 

W ( i , 1 ) =a ( i , 1 ) *w ( i , 2 ) +a ( i , 2 ) *w ( i , 3 ) +e ( j ) 
xbp ( j ) =b ( i , 1 ) *w ( i , l ) +b ( i , 2 ) *w ( i , 3 ) 
e( j+l)=xbp(j)+e(j) 
i=i-l 

10 continue 
en=e (n+1) 
c 
c 

i=n 

do 20 k=l , n 
u (k) =en-xbp (k) 

g(i,3)=g(i,2) 

g(i,2)=g(i,l) 

g(i,l)=a(i,l)*g(i f 2)+a(i,2)*g(i,3)+u(k) 
gbp (k) =b ( i , 1) *g (i , 1) +b ( i , 2) *g ( i , 3) 
q(k, 3)=q(k, 2) 
q(k, 2) =q(k, 1) 

q(k, l)=a(i, 1) *q(k,2)+a(i,2) *q(k,3) +gbp (k) 
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s(i)=c(i)*q(k,2) 

i=i-l 

20 continue 

kk=mod (ki , 10) 

c if(kk.eq.O) write(9,*) ki, en, a (1 , 1) , as (l) , a (2 , 1) , as (2 ) 
c if(kk.eq.O) print *, ki, (a (i, 1) , i=l,n) 

c if(kk.eq.O) print *, ki, (as (i) , i=l , n) 

c if (ki.ge. 1800) write(9,*) ki , e (1) , en,xbp(l) , a (1 , 1) , as (1) 
c if (ki .ge. 1800) print *, ki, e (1) , en, xbp (1) , a ( 1, 1) , as (1) 

return 
end 
c 
c 
c 

subroutine system (yk, wk, theta, step, var, ic, n, k) 
dimension theta(10) ,yout(10) ,yz(10) 
if (k.ne.l) goto 100 
100 continue 

call bpf (theta , yz , n, k) 
call waves (theta, yout,n,k) 
call file (yk, theta, step, ic, n, k) 
sum=0. 0 
wk=0 . 0 
do 10 i=l,n 
sum=sum+yout ( i ) 
wk=wk+yz (i) 

10 continue 

call gnoise (g, var, k) 
c yk=sum+g 

c yk=wk+g 

yk=yk+g 
return 
end 
c 
c 

subroutine waves (theta, yout,n,k) 

dimension zkar(10,3) ,theta(10) ,yout(10) ,al(10) 

if (k.ne.l) goto 10 

pi=(atan(l. 0) ) *4 . 0 

rad=pi/180 . 0 

do 20 i=l,n 

zkar(i, 1)=0.0 

zkar ( i, 2 ) =-sin( theta (i) *rad) 
al (i) =2 . 0*cos (theta ( i) *rad) 

20 continue 
10 continue 

do 30 i=l,n 

zkar (i , 3 ) =zkar ( i , 2) 

zkar (i , 2) =zkar (i , 1) 

zkar(i, l)=al(i) *zkar(i,2) -zkar(i, 3) 
yout ( i) =zkar ( i , 1) 

30 continue 
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return 

end 

c 

c 

subroutine gnoise (g, s, k) 
if(k.ne.l) goto 10 
ix=l 

yfl=rand ( ix) 
ix=0 

10 continue 
sum=0 . 0 
do 20 i=l, 12 
yfl=rand(ix) 

20 sum=sum+yfl 

g= (sum-6 . 0) *s 
return 
end 
c 
c 

subroutine stability (akl , ak2 , flag, n) 
dimension akl (10) ,ak2 (10) , flag (10) ,a(10,2) 
do 30 i=l,n 
flag(i) =0. 0 

a (i , 1) =- (akl (i) *akl (i) +ak2 (i) -2 . 0) 
a ( i , 2 ) =- ( 1 . 0-ak2 ( i) ) 
des=a (i , 1) *a (i , 1) +4 . 0*a (i , 2) 
if (des.lt. 0.0) goto 10 
des=0 . 5*sqrt (des) 
rootl=0 . 5*a ( i , 1) +des 
rl=abs (rootl) 
root 2=0. 5*a (i, 1) -des 
r2=abs (root2) 
if (rl.ge.1.0) goto 20 
if (r2.ge.l.O) goto 20 
c print *, '-real-roots — >', rootl , root2 
goto 30 

20 flag(i)=1.0 

goto 30 
10 des=-des 

y=0 . 5*sqrt (des) 
x=0 . 5*a ( i , 1) 
xx=abs (x) 

if (xx. le . 1 . Oe-6) x=0. 0000001 
angle=(atan(y/x) ) *57 . 3 
radius=sqrt (x*x+y*y ) 

c print *, ' — complex — roots — >', radius, angle 
if (radius. ge. 1. 0) goto 20 
30 continue 
return 
end 
c 
c 



81 



c $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 

c 

c 

subroutine bpf (theta, yout,n,k) 

dimension zkar(10,3) ,a(10,3) ,zkma(10,4) ,b(10,3) 
dimension theta(lO) ,g(10) ,yin(10) ,yout(10) ,uk(10) 
if (k.ne.l) go to 10 
pi=4 . 0*atan (1.0) 
rad=pi/180 . 0 
do 15 i=l,n 
zkar (i , 1) =0 . 0 
15 zkar (i, 2) =0. 0 
r=0 .985 

c print *,' — input — r — >',r 

c read *,r 

do 25 i=l,n 

a (i, 1) =1 

a(i,3)=r*r 

g ( i) =0 . 5* ( l-r*r) 

(i , 1) =g ( i) 

b(i,3)=-g(i) 

a (i , 2 ) =2*r*cos (theta (i) *rad) 

25 b ( i , 2 ) =0 . 0 

10 continue 
save=var 
do 35 i=l,n 
var=2 . 0 

call gnoise (gk,var ,k) 
yin ( i) =gk 

zkar ( i , 3 ) =zkar ( i , 2 ) 
zkar ( i , 2 ) =zkar ( i , 1) 
zkma ( i , 3 ) =zkma ( i , 2 ) 
zkma ( i , 2) =zkma (i , 1) 
zkma ( i , 1) =yin (i) 

uk(i) =zkma ( i , 1) *b( i , 1) -zkma ( i , 2) *b (i , 2) +zkma (i, 3 ) *b ( i , 3) 
zkar (i, 1) =zkar (i, 2) *a (i, 2) -zkar (i, 3) *a (i, 3) +uk(i) 
c print ' — fil — >' ,afil,theta,r 

c print — fil — >',zkar 

35 yout ( i) =zkar (i , 1) 

var=save 
return 
end 

c $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 

c 

subroutine f ile (yk, theta, step, ic, n, k) 
dimension output (8000) , theta ( 10) 
if (k.ne.l) goto 10 

call gendata ( output , theta , step , ic , n , k) 

10 continue 

yk=output (k) 

return 

end 
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c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



c 

subroutine gendata ( output , theta , step , ic , nn , k) 

DIMENSION YSAMP ( 65536) , data (4 , 8000) , output (8000) 
dimension theta (10) 

c OPEN (UNIT=9 , FILE= ' bpsk.dat ' , STATUS= 'NEW ' ) 
n=7000 
amag=l .414 
tcarr=4 . 8 
tchip=7 . 0 
ichip=0 
tdata=n 
tdelay=0 . 0 
1=9999 
tchip=200 . 0 
do 10 i=l,nn 
tcarr=360 . 0/theta ( i) 

call bpsk (n, amag, tcarr , tchip, ichip, tdata , tdelay, 1 , ysamp) 

do 20 j=l,n 

if (ic.eq.0) goto 80 

if (i.eq.2) goto 50 

data ( i , j ) =ysamp ( j ) 

goto 60 

50 if (j.gt.3000) ysamp(j)=0.0 
80 data ( i, j ) =ysamp ( j ) 

60 continue 
20 continue 

print *,i,' — >', theta ( i) , tcarr 
10 continue 

tcarr=360 . 0/90.0 
tchip=8 . 0 

call bpsk (n , amag, tcarr , tchip , ichip , tdata , tdelay, 1 , ysamp) 
do 30 i=l,n 
x=0 . 0 

do 40 j=l,nn 
40 x=x+data ( j , i) 

output ( i ) =x+y samp ( i ) 
c write (9,*) x 
30 continue 

if (ic.eq.0) return 
nl=n-3000 
tchip=200 . 0 
theta l=theta ( 2 ) +step 
tcarr=3 60 . 0/thetal 

call bpsk ( nl , amag , tcarr , tchip , ichip , tdata , tdelay , 1 , ysamp) 

j=l 

do 70 i=3000, n 
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output ( i ) =output ( i ) +ysamp ( j ) 
j=j + l 

70 continue 

c close(9) 

return 
end 
c 
c 

c this programme was modified to suit the martin 

c programme 

subroutine bpsk 

& (n, amag, tcarr , tchip, ichip, tdata , tdelay , 11 ,ysamp) 



c 

c 

C THIS PROGRAM GENERATES SAMPLES OF A DIRECT SEQUENCE BI-PHASE 
SHIFT 

C KEYED SPREAD SPECTRUM SIGNAL. THE "INFORMATION" BITS AND THE 
C SPREADING SEQUENCE USED ARE RANDOMLY GENERATED. PARAMETERS 
REQUIRED 

C FOR OPERATION ARE: 

C N = NUMBER OF SAMPLES GENERATED 

C FOR CONSISTENCY WITH USE BY AN FFT 

C ALGORITHM, N SHOULD BE AN INTEGER 

C POWER OF 2 — TYPICALLY 1024 

C NOTE: IF N>1024 DIM OF YSAMP MUST CHANGE 

C MAG = MAGNITUDE OF CARRIER WAVEFORM 

C TSDELAY = NUMBER OF SAMPLE TIMES DELAYED FROM 

C t = 0 BEFORE BEGINNING SAMPLES 

C AN ARBITRARY VALUE THAT ALLOWS SOME 

C FLEXIBILITY IN SAMPLING. SHOULD BE 

C ZERO FOR SAMPLING AT t=0. 

C TDATA = NUMBER OF SAMPLE TIMES IN ONE DATA BIT 

C TCHIP = NUMBER OF SAMPLE TIMES IN ONE BIT OF 

C SPREADING CODE 

C TCARR = NUMBER OF SAMPLE TIMES IN ONE CYCLE OF 

C CARRIER FREQUENCY 

C TSAMP = DURATION OF A SAMPLE TIME. IN GENERAL 

C TSAMP WILL ALWAYS BE = 1.0, SINCE 

C TIME SAMPLED VALUES BECOME STATIC 

C WHEN STORED AND CAN BE SCALED LATER. 

C 

C THE ALGORITHM USED TO GENERATE THE DS-BPSK SIGNAL RELIES UPON 
THE 

C BUILTIN FORTRAN RANDOM NUMBER GENERATOR RAND(L) TO PRODUCE THE 



C INFORMATION BIT STREAM AND THE SPREADING SEQUENCE CODE. AT 
THE TIME 

C OF EACH SAMPLE, THESE TWO BITS ARE MODULO 2 ADDED TO PRODUCE 
A CHIP 

C BIT. THIS BIT IS THEN USED TO SHIFT THE PHASE OF THE CARRIER 
BY 
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C 180 DEGREES IF 1 OR BY 0 DEGREES IF 0. THE RESULTANT VALUE OF 
THE 

C COSINE FUNCTION IS MULTIPLIED BY THE AMPLITUDE OF THE 
WAVEFORM, THEN 

C STORED IN FILE BT.DAT IN THIS FILE DIRECTORY. THE PROCESS IS 
CONTIN- 

C UED UNTIL N SAMPLES ARE GENERATED. NOTE THAT FIRST LINE OF 
BT.DAT 

C CONTAINS N, DSBPSK, DATE AND TIME AS ID FIELD FOR FILE. 

C 

C W.R. TUCKER 9-MAY-83 

C CONVERTED PROGRAM TO ALLOW FOR M-SEQUENCE CHIP 

C M-SEQUENCE GENERATOR REQUIRES DATA FILE FSR.DAT WITH 

C PARAMETERS FOR FEEDBACK SHIFT REGISTER. 

C W.R. TUCKER 4 -AUG-8 3 

C 

C INITIALIZE — I = INFO BIT NUMBER, J = SP SEQ CODE BIT NUMBER 

C L IS ARGUMENT FOR RAND(L) — A DIFFERENT SEQUENCE MAY BE 

C GENERATED BY INITIALIZING L WITH DIFFERENT VALUES. 

C 

C Modified by HHL for installation on ULTRIX ECE VAX, 12/90 
c Writes output in file bpsk.dat 
c 

DIMENSION YSAMP (65536) 

INTEGER I , J , L, INFO , CHIP, IFSR, IFBD (100) ,IVAL(100) 

REAL TDATA, TCHIP, TCARR, TDELAY , TS AMP, MAG , YARG , YSAMP, PI 
PI = 3.14159265 
TSAMP =1.0 
1 = 0 
J = 0 
L = 0 

WRITE ( 6 , 900) 

900 FORMAT ( ' GENERATION OF DIRECT SEQUENCE BPSK SIGNAL', 

A ' IN FILE bpsk.dat') 

WRITE (6, 1000) 

1000 FORMAT ( ' # SAMPLES TO GENERATE = ',$) 

c READ ( 5 , * ) N 

print *,n 
WRITE (6, 1001) 

1001 FORMAT (' MAX AMPLITUDE OF SAMPLE VALUES (R)= ',$) 

c READ ( 5 , *) MAG 

mag=amag 
print *,mag 
WRITE (6, 1002) 

1002 FORMAT (' # SAMPLES PER CARRIER CYCLE (R)= ',$) 

C READ (5, *) TCARR 

print *,tcarr 
WRITE (6, 1003) 

1003 FORMAT (' # SAMPLES PER CHIP BIT (R) = ',$) 

C READ (5, *) TCHIP 

print *,tchip 
WRITE(6, 1020) 
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1020 FORMAT ( ' ENTER ( 0) : RANDOM CHIP, OR ( 1) : REPEAT M-SEQ CHIP 
= ' / $ ) 

c READ ( 5 , *) ICHIP 

print *,ichip 
WRITE (6, 1004) 

1004 FORMAT (' # SAMPLES PER INFO BIT (R) = ',$) 

c READ (5, *) TDATA 

print *,tdata 
WRITE(6, 1005) 

1005 FORMAT (' # DELAYS BEFORE SAMPLING (R)= ',$) 

c READ (5, * ) TDELAY 

print *,tdelay 
WRITE(6, 1006) 

1006 FORMAT (' RANDOM NUMBER SEED (14) > ',$) 

1=11 

C READ (5, * ) L 

print *,1 

c Initialize RAND by calling with input non-zero L. 

c Subsequent calls will be with L = 0. 

RANDOM=RAND ( L) 
c debug 

C WRITE (6,*) ' SEED AND RANDOM NUMBER RETURNED L, RANDOM 
L=0 

WRITE(6, 1010) 

1010 FORMAT (/IX, ' WORKING ' ) 

IF (ICHIP .EQ. 0) GO TO 5 

OPEN (UNIT = 1, FILE = ' FSR. DAT ' , STATUS = 'OLD') 

READ (1,500) IFSR 

READ (1,600) (IFBD(K) , K = 1, IFSR) 

500 FORMAT (13) 

600 FORMAT (13) 

CLOSE (UNIT = 1) 

DO 5 K = 1 , IFSR 
IVAL(K) = 1 
5 CONTINUE 

DO 100 K = 1 , N 

C CHECK TO SEE IF WE NEED TO GENERATE A NEW DATA BIT 
IF ( (K+TDELAY) *TSAMP .LT. I *TDATA*TSAMP) GO TO 10 
C IF SO DO IT HERE 

RANDOM = RAND(L) 
c debug 

C WRITE ( 6 , * ) • IFLAG AND RANDOM NUMBER RETURNED L, RANDOM 

INFO = 0 

IF (RANDOM .GT. 0.5) INFO = 1 
C KEEP TRACK OF WHICH DATA BIT WE ARE ON 
I = 1 + 1 
10 CONTINUE 

C NOW CHECK TO SEE IF WE NEED TO GENERATE A NEW CHIP BIT 
IF ( (K+TDELAY) *TSAMP .LT. J*TCHIP*TSAMP) GO TO 20 
IF (ICHIP .NE. 1) GO TO 15 
CHIP = IVAL(IFSR) 

CALL MSEQ ( IFSR , IFBD , IVAL) 
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GO TO 19 
15 CONTINUE 
C IF SO DO IT HERE 

RANDOM = RAND (L) 

CHIP = 0 

IF (RANDOM .GT. 0.5) CHIP = 1 
C KEEP TRACK OF THE CHIP BIT NUMBER 

19 CONTINUE 
J = J+l 

20 CONTINUE 

C NOW WE DETERMINE THE PHASE SHIFT 

PHASE = FLOAT (MOD (INFO + CHIP, 2)) 

C COMPUTE THE ARGUMENT FOR THE COSINE FUNCTION 

YARG = 2 . 0*PI* (1.0/ (TCARR*TSAMP) ) * (K+TDELAY) + PI*PHASE 
YSAMP(K) = MAG * COS (YARG) 

100 CONTINUE 
C NOW SAVE THE VALUE 

c ranga OPEN (UNIT=1,FILE=' bpsk.dat ' , STATUS= 'NEW ' ) 

C CALL DATE ( BDATE ) 

c CALL TIME (CTIME) 

c-ranga WRITE ( 1 , 12 5 ) N , MAG , TCARR , TCHIP , TDATA 

125 FORMAT ( 15 , 5X , F4 . 1 , 10H*DSBPSK+M2 , F7 . 2 , lHf , F7 . 2 , 1HC , F7 . 1 , 

* lHd) 

c-ranga WRITE(1,150) (YSAMP(I) , 1=1 ,N) 

do 121 i=l,n 
c write (1,*) ysamp(i) 

121 continue 

150 FORMAT (8F16. 8) 

CLOSE (UNIT=1 ) 

return 

END 

SUBROUTINE MSEQ (IFSR, IFDB, IVAL) 

C THIS SUBROUTINE PERFORMS THE SHIFTING OPERATION OF AN IFSR 
STAGE 

C FEEDBACK SHIFT REGISTER, WITH FEED BACK CONNECTIONS AS 
INDICATED 

C BY IFDB. IVAL IS THE INITIAL CONTENTS OF THE FSR AND WILL 
CONTAIN 

C THE FINAL CONTENTS AFTER SHIFTING. 

C W.R. TUCKER 4 AUG 83 

INTEGER IFDB (IFSR) , IVAL(IFSR) 

I SUM = 0 

DO 10 I = 1, IFSR 
ISUM = IFDB (I) * IVAL ( I ) + ISUM 
10 CONTINUE 

IBIT = MOD ( ISUM, 2 ) 

DO 20 I = 1, IFSR - 1 

IVAL (IFSR + 1 - I) = IVAL (IFSR - I ) 

20 CONTINUE 
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IVAL(l) = IBIT 

RETURN 

END 

c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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